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ABSTRACT 


An algorithm for the identification of the input distribution 
matrix of a linear, stationary system operating in a stochastic 
environment is derived. The identification is accomplished by defining 
a set of autocorrelation functions for a noise element composed of a 
linear combination of the input distribution matrix elements and the 
random excitations of the system. Another possible identification 
method employing a Kalman filter is discussed and the problems associ- 
ated with its derivation are presented. Results of computer simulations 


for both methods are included. 
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I. INTRODUCTION 


In recent years, an emphasis on plant identification for a system 
operating in a stochastic environment has developed in the field of 
automatic control. It is necessary to know the dynamics of the system 
before it may be effectively controlled. In order to know the system 
dynamics, a suitable mathematical model of the system must be found. 
For a system subjected to a random forcing function, finding a suitable 
mathematical model involves the identification of both the state 
transition matrix and the input distribution matrix for the random 
forcing function. Several methods have been proposed for the identi- 
fication of the state transition matrix, but no one has proposed a 
method for identifying the input distribution matrix for systems with 
a single random input. The purpose of this investigation is to propose 
such a method. The problem may be defined as follows: 

Given: 

1. <A system whose behavior may be defined by a set of linear, 
constant-coefficient, differential equations. 

2. A statistical description of the excitation of the system. 

3. A sequence of observations on the state vector. 

4. An estimate of the state transition matrix of the system. 

Problem: 

Determine the input distribution matrix of the system. 

Although this investigation is concerned only with systems with 
random forcing functions, the methods described may be extended to 


Systems with both random and control inputs. The method of section 
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III may be extended by changing the form of the estimated state transi- 
tion matrix as described by lee [5]; and that of section IV, by including 
the effect of the control input in the filter equations as described by 
Kalman [3]. The problem is restricted to random forcing functions for 


simplicity. 


WZ 


IIT. BASIC MATHEMATICAL MODEL 


The investigation reported here starts with the assumption that 
there exists a set of n linear, constant-coefficient, differential 


equations of the form 


he 


x(t) = Ax(t) + Bw(t) 1) 
which describe the dynamic behavior of a system subjected to a random 


forcing function, w(t). The solution of these equations is 
t 


x(t) = ?(t-to) x(0) *f (t-T) Bw(t)dtT (222) 


t 
where 0 


B(t) =f"? ((st - ay. 


Because the continuous system dynamics will be observed at discrete 
times, it is necessary to describe the system by discrete or difference 
equations. Introducing a sampling device of period T and a zero- 
order Weld on the forcing function allows the system equations to be 


written 


eS)er4] = 2 (T) xy + U(T)wy (2.2) 


where x; is the system output vector at a time kT, ® is the nxn discrete 
State transition matrix, | (T) is the 1lxn input distribution matrix and 
wy is the sampled and zero-order held input vector at time kT. The 


components of the input vector w, are assumed to form a gaussian white 


k 
sequence of zero mean and known variance during the period of system 
operation under investigation. 

Observations made on the system output are assumed to be scalar 


linear combinations of the states of the system which are contaminated 


by additive gaussian white noise of zero mean and known variance. At 


th . ; 
the k sampling instant the observation z, is given as 


k 
a = Hix, ate VE (2.4) 

where H is the lxn known measurement matrix and vie is the scalar 
additive measurement noise. The system may then be represented by 
the block diagram of Figure l. 

~ In the general problem of identification of linear dynamic systems, 
it is necessary to estimate the elements of the state transition matrix 
before proceeding to identify the input distribution matrix. For the 
purposes of this investigation, this estimation is assumed done by 
using the method described by Lee [5]. The estimation of the state 
transition matrix by this method makes it necessary to change the 
system model to the so-called canonic form. Assuming that the system 


is observable and n-identifiable , the difference equation for the 


system model then becomes 
= Ox* * 
eae. eae 


(223) 


where 


(236) 
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Dk¥ = Cl = b = 2 


H* = { ii O O eee O ] 
The matrix is now the estimated state transition matrix obtained by 
using the method of Lee. The canonic system model can then be 


represented by the block diagram of Figure 2. 





ve 3 a ; k+1 “ae =k me) 
d 


FIGURE 1. BLOCK DIAGRAM OF THE SYSTEM 








FIGURE 2. BLOCK DIAGRAM OF THE 
SYSTEM IN CANONIC FORM 
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III. DERIVATION OF THE IDENTIFICATION ALGORITHM 


Consider a system with canonic state equations of the form 


J) 5 ees 
Coe) 
ee gel 
where Wy is assumed to be a scalar random forcing function. The 
z-transformed state equation becomes 
zy(z) - zy(0) = 6* y(z) + P* wz). (3.2) 


Assuming that the initial conditions y(0) are zero or that the 
identification process does not start until the initial conditions 
no longer affect the system output, the system may be represented by 
the signal-flow graph of Figure 3. From this signal-flow graph, the 
transfer function Z\z)/W(z) may be written as 

b b oe | uy ? 

Bega -2)* a Oe ee 

Z(z) _Z Z Z Zz 

ret 2 a 








Dae 
n-2 


(3339) 





Since the elements of the ®* matrix are assumed to be known, the 
problem can now be stated as finding the location of the zeros of 


the transfer function and thus the values of bi» bo > ee bo 


From the transfer function, the difference equation relating the 


measured output z, to the input w, may be written as 


k k 
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i Pian Potent’) Once oes ee 


eS eee (3.4) 

5 OE ei eee 

- Ss ie i oP eee 
Shifting indices by +1 and using vector notation, Equation (3.4) 
becomes Za] Pa Cade nid ac (569) 
where 





Gils) 





ee 


Because the elements of both the vectors n 


ie and d are unknown, a single 


noise element is defined as 


ndtv, (3.7) 


where Oy. 1s a sca leterand 


tl 
Ou. 


O 1x * 4 Taig oo Cam 


Toes OT re 7 ec) ea 


(35a) 


ects act 


! 
a. 


ra | ae ne ee Peete oe 


From Equation (3.8) it is apparent that the noise element oF is 
correlated with the elements OnnE OD? ree hen T3 however, oF 
and Qn are independent with respect to each other. This property 
of the noise elements is the basis for the identification algorithm. 
As shown in Appendix A, a set of autocorrelation functions for the 


noise elements may be defined such that 


2 2 2 2 2 2 
(0) 0, Roe GL cl aP oa ae lm) = 9 i 
(3.9) 
ae 
RG@py— 6 [9,04 7! a0 a Hes j=l, 2, 2.2mm 
R(j) = 0 j>n 


where ome is the standard deviation of the gaussian white sequence of 
noise inputs Wy and om is the standard deviation of the measurement 


noise. 


By rearranging Equation (3.5) such that 


7 T 
Same carer) Sy, O , (3.10) 
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the noise element is now a linear function of the system outputs 

a Zn? ae 2] and is therefore known. Hence, the auto- 

correlation functions of equation (3.9) may be found by taking the 
any: 2 

Miatistics £ [Q. Jere LOO g Equation (3.9) now represents a set 


of n simultaneous nonlinear algebraic equations which may be solved 


ane values of d,> d,» Nees di: The values of b> bo» Sites DO may 


then be found from equation (3.6). 

The general procedure for the identification of the input distri- 
bution matrix in a linear system with a random forcing function may be 
outlined as follows: 

1. Measure a scalar linear combination of the states of the system 
plus additive noise. 

ip Identify the canonic form of the state transition matrix of 
the system. | 

3. Galculate the values of the noise element oF at each sampling 
time. 

wee catculate the autocorrelation functions R(j) for j] = Oe...) 
n-l. 

5. Solve the simultaneous equations (3.9) for the elements of d. 

6. Solve equation (3.6) for the elements of b. 

Although the above algorithm does provide a means for identifying 
the input distribution matrix, two problems arise in the identification 
process. The first of these involves deciding when to start the 
identification process. Because the autocorrelation functions are 
computed by taking the average of the values of - Oa 7» —, 

OO anHl at a number of sampling instants, this process must not begin 


until enough time has elapsed so that the sample mean and variance of 


ZL 


; ; Z 
the noise input and the measurement noise approach 0, oO and 0, o 
W V 


respectively. If the identification process is begun too soon, the 
computed values of b will be in error from the actual values. 

The other problem in the identification process can be best shown 
by an example. Assuming n = 2, the autocorrelation functions then 
become 

Z 


om (d, 


Z 


it 


Z Z 
R(O) + d, y+ on 


(3S 


2 
R(1l) = oO. dy d,, : 


Solution of Equation (3.11) gives 








20 
Ww 
(3. ty 
eZ 1/2 
OW 
and from Equation (3.6) 
b, = 4, 
(3. 135 


Oy yO eas 


It is apparent from Equation (3.12) that the identification of the 
elements Di: bo involves a decision as to which combination of signs 
to use as the correct one. For higher order systems this decision 


becomes more complex since the number of possible sign combinations 


increases. 


ce 


IV. KALMAN FILTER APPROACH 


An alternative to the identification algorithm of the preceding 
section may be derived from the equations of the Kalman filter 
described below. Consider again a system with the canonic state 
equations of the form 


Let] k 


(4.1) 


li 


we 
Zc] 1 eS 


where w, is a scalar random forcing function. The canonic state 


k 
vector y, may be estimated by the technique described by Kalman [3] 
where, given k previous measurements, the estimate vine at the ae 
samp ling instant is given by 
= met sa - H*® 6% § 
Dek Preeti aoe : ye (4.2) 


where G. is the filter weighting or gain applied at the eee iteration 


and H* is the lxn known measurement matrix. The gain is calculated 


from a set of recursive equations of the form 


i T ~1 
= 4e one 
G. Ween H (H* PL fk=1 He + R] 
= - * 
Ei tS E") Een (4.3) 
p =$*p, 4° 4Q 
k+1/k k/k 
where Pk and Poti /k are the estimation and prediction covariances, 
respectively, and R and Q are given by 
A 
R=E [vy] 
A 9 7 (4.4) 
QTR E [wo] Te, 


fae 


Since all elements necessary to predict the canonic state vector 
y, are known except the input distribution matrix [*, the system model 
based on this estimation scheme will produce suitable estimates for By 
only if the correct [* is used in the filter equations. In order to 
determine if the correct !'* is being used, it is necessary to define a 
performance index to compare the performance of the model with that of 


the system. Because the actual system output is available only in the 


form of the measurement Zo the performance index must be defined in 
' N : 
terms of Ze and the predicted measurement Zr N=1 as given by 
A A 
= Fe ox 
jc RN oT oy 


which is obtained from the system model. 

Before defining the performance index to be used, it is necessary 
to define some other terms. The filter is said to have reached the 
steady-state condition when the prediction covariance matrix does not 
change from one sampling instant to the next. Under the steady-state 
condition the actual prediction measurement variance is given by 
K+M re 9 


A ey ee Cie 
i=K 


EZ(M) = 


<a 


where K is the value of k at the sampling instant when the filter 
reaches the steady-state condition. The integer M is chosen such 
that the relation 

|Ez(M+1) - EZ(M)| < 0.05|Ez(™)| (4.7) 


is valid. The predicted measurement variance is defined as 


A 2 


pz = Sine ey ee 


(4.8) 
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which, as shown in Appendix B, may be written as 


PZ = He PL qx +R (4.9) 


where Pe is the prediction covariance under the steady-state conditions 
of the filter, 

From the above definitions the performance index J is defined as 
the difference between the actual and the predicted measurement 


variances and may be written as 


J = EZ(M) - PZ 


1 K+M 


ae Ee o iia ns ; c *K He + | ° (4.10) 


> 


wemce it is desired that the output of the model 2 / and the output 


k/k-1 


of the system z, be equal for perfect operation of the model, J is to 


k 
be minimized. Because both EZ(M) and PZ are functions of the input 
aistribution matrix I*, the minimum for J will result when the correet 
['* is used in the model. The problem may now be stated for this 
approach as identifying I* so as to minimize the performance index J. 
Because it would waste time and effort, a trial-and-error search 
femeethe correct |[* is undesirable. Therefore, it is necessary to 
derive some method of changing the elements of I'* used in the model 
after ech calculation of J, based on the value of J. The method 
investigated is based on a gradient search where the elements of ['* 
are changed by some amount according to the values of a sensitivity 
function and J. The sensitivity function must reflect the change in 


J caused by changing I[*. The method of changing I’* may be written as 


3c = | *¥ 
re = Te + £ CS, TH) J (4.11) 


25 


where the subscripts "old" and "new" represent the value of I’* just 
used in the model and the value to be tried next, respectively. The 
sensitivity function is denoted by f£ (J, ['*) where f£ is a vector 
Pune t Lon. 

The problem is now to find a suitable sensitivity function to use. 
Two functions and their failures during simulation are presented in 
the next section. Because the performance index is a function of the 


measurements 2 (which are functions of all past 


K? KHL? °°? 7KM 
measurements) and the predicted measurements during the interval K, 
K+M, (which are also functions of all past measurements), the problem 


of finding a suitable sensitivity function is a complicated one. At 


present no suitable function has been found. 
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V. COMPUTER SIMULATIONS 


In order to determine the accuracy of the proposed identification 
algorithm, a second-order system was simulated on the IBM 360/67 com- 
puter. For the system chosen the canonic state transition matrix, the 


input distribution matrix, and the measurement matrix were 


0 1 
6x = 
=f 16. lie 
i 
[* = (st) 
U7 2 
H* = [1 0 | ; 


In the simulation, the actual and the canonic-state transition matrices 
were assumed equal so that no error in the identification process would 
result from a difference in the two,and also for simplicity. The random 
oneine function and the measurement noise were obtained from a random 
number generator which gives a gaussian white sequence of any desired 
mean and standard deviation. The means and standard deviations for the 


two sequences were 


E [w = 0 
Z 
E [w, J] = 2 
(5.2) 
E [v,] = Q 
E v,") = le 


The system dynamics were simulated and the measurements Zz were used 


to obtain the autocorrelation functions R(O) and R(1). These values 
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were "then used to compute the values of by and bo at each sampling 


instant from the equations 


2 ie 2 
b, (k) = | 0.125 | RO) Sel ta ARO a= be = 4R(1)*] j | 


(5.85) 


1/2 51/2 
by (k) = 0.125 | R(0) - 1 - [(R(O) - bi - 4R(1)°) in 


The system and the identification process were simulated for 1000 
sampling intervals of period one second. An ensemble of 200 simu- 
lations was run and an ensemble average for R(O), R(1), b, kk), and 

by (k) was taken. The results at various sampling intervals are shown 
in Table I. The ensemble averages of the differences between the 
actual values of I’* and the computed values are also shown. 

As can be seen from Table I, there is a bias in the autocorrelation 
function R(1) which produces a bias in the estimates for by and bo: 
No explanation for the bias can be given. Simulation of another, 
Similar second-order system also produced a bias, equal to -1.0. 

For the algorithm to produce the desired results, it is necessary to 
eliminate the bias. 

A second-order system was also simulated for the Kalman-filter 


approach to the identification problem. The values used for the 


simulation were 


0 1 
Ok = 
1/8 1/ 
1 
T% = (5.4) 
1/2 
H* = [1 0 | 
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The random forcing function and the measurement noise were generated as 


described above with 


E [w,] = 0 
E fw] = 3 
E [v,] = 0 
B[v, ]=R=4 « 


(5.5) 


The performance index J was calculated for various values of I* with 


the results given in Table II1. It can be seen that J is a minimum 


when the correct I* is used. 


Two sensitivity functions were also calculated during the simulation. 


The first as derived in Appendix C is given by 


_ be 
£) 2 2%) Spee 


The second, also derived in Appendix C,is given by 


A )? | ByoEZ 


S ws : OP4 
£, Gi) = WY orp) = Gem 25s or 


As shown in Table II, neither of these sensitivity functions is 


suitable since neither produces a consistent pattern by which to 


(Ec) 


(567) 


change T*, The failure of these functions appears to lie in the fact 


that neither measures the sensitivity of the measurements to the values 


ee |, 
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: Z 
k b, (k) bby (k) by (k) 2 by (kK) BIG Ie ELA M4] 


2 0.0 0.0 0.0 0.0 5.11 0.0 

10 0.84 0.16 0.65 =f), 15 5.25 -0.31 
100 1.003 -0.003 0.645 aM Bild -0.49 
200 1.009 -0.009 0.635 0.135 5.18 -0.50 
300 1.010 -0.01 0.633 138 iy -0.50 
400 1.012 -0.013 0.633 00133 5.18 -0.50 
500 1.016 -0.016 age 20.131 5.20 -0.50 
600 1.017 -0.017 0.634 -0.134 5.21 -0.51 
700 1.016 -0.016 0.633 miigK 5.20 Si 
800 1.018 -0.018 p34 -0.135 5.22 One 
900 1.019 -0.019 0.636 -0.136 5.23 Bi 
999 1.019 -0.019 0.635 70.135 5.23 Bi) 
RCO) = oman + d.°) + 0. = 5.0 (should match E[0,“1) 
R(1) = 0,” dd, = 0.0 (should match E[9,0,,,1) 

W iL k kt1 


TABLE I. Results of Simulation 
of Identification Algorithm 
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*Correct I* 


EZ (M) 





14.7 
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14.3 
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oie 
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Kalman Filter Approach 
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Results of Simulation of 
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4, -0.6 
ate 2 
i, 20n3 
1, 0.04 
5, -0.8 
7, aie 
ES 
0,-28 lo 
I ote 


VI. CONCLUSIONS 


The identification of the input distribution matrix of a linear 
system operating in a stochastic environment is a complicated process, 
but it can be done to some degree by the proposed algorithm. Further 
investigation is needed to determine the reason for the bias in the 
autocorrelation function R(1). A general procedure for choosing the 
proper signs to use in the solution of the nonlinear algebraic equations 
is also desirable. 

The Kalman-filter approach would be an excellent method to use in 
the identification process if a suitable sensitivity function can be 
found. This function must reflect the sensitivity of the performance 
index to the changes in the input distribution matrix. A possible 
first step in finding the function is to examine the sensitivity of 


themmeasubrements towthe elements Of ginesmatrix. 
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DERIVATION OF THE AUTOCORRELATION FUNCTIONS 


APPENDIX A 


FOR THE NOISE ELEMENT 
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APPENDIX B 


DERIVATION OF THE PREDICTED MEASUREMENT VARIANCE 
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APPENDIX C 


DERIVATIONS OF THE SENSITIVITY FUNCTIONS 
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The first sensitivity function to be derived is an during steady- 
state operation of the filter. This is done by finding an expression 
for the prediction covariance during steady-state operation, equating 
elements on both sides of the equation, and then taking the desired 


partial derivatives of the implicit functions. The derivation follows. 
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Another sensitivity function of the form 
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is to be considered now. It may be found by expressing the filter gain 
in steady state as a function of the steady-state prediction covariance. 


The required partial derivatives may then be found as above. The follow- 


ing derivation is done assuming steady-state operation of the filter. 
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where re is as given above. 
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